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Following the work of Leutheusser [Physica A 127, 667 (1984)], the solution to the Percus- Yevick 
equation for a seven-dimensional hard-sphere fluid is explicitly found. This allows the derivation 
of the equation of state for the fluid taking both the virial and the compressibility routes. An 
analysis of the virial coefficients and the determination of the radius of convergence of the virial 
series are carried out. Molecular dynamics simulations of the same system are also performed and a 
comparison between the simulation results for the compressibility factor and theoretical expressions 
for the same quantity is presented. 

I. INTRODUCTION 

In liquid theory there has been a long lasting interest on the equilibrium properties of high-dimensional hard-sphere 
fluids, especially in the last few years. ^^'^'^ Such an interest has arisen from many different sources. To begin with, 
given the relative simplicity of the intermolecular interactions in these hard-core systems, they are amenable to both 
theoretical and computer simulation studies. In this sense and as it occurs in other problems in theoretical and 
mathematical physics, it is an asset that one can deal with hard spheres in arbitrary dimensionality and exploit some 
of the features that these systems have in common, for instance the fact that they all exhibit a first-order freezing 
transition. Furthermore, and as conjectured by Frisch and Percus^ in the case of fluids high spatial dimensionality 
may have a parallel with limiting high density situations, so that by increasing the dimensionality one may obtain 
at least a rough idea of any thermodynamic phenomenology that extends to such dimensionality. An example of this 
expectation is the recent investigation of the demixing problem in mixtures of hard hypersphereSf?- 

Computer simulation studies of hard-sphere fluids in dimensions greater than three are very scarce. To the best 
of our knowledge, only the four- and five-dimensional simple-'^- and multicomponent^^ fiuids have been simulated. 
This is not surprising since the computational effort needed to obtain reliable results increases significantly with the 
dimensionality. 

Exact information on the equation of state (EOS) usually comes from the virial coefficients Bn defined bji^i 
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In this equation, p is the pressure, p is the number density, fcs is the Boltzmann constant, T is the temperature, and 
Z is the compressibility factor. The second virial coefficient is B2 — 2'^~^t;d<T'^, where d is the dimensionality, a is 
the diameter of a sphere, and Vd = (7r/4)''/^/r(l + d/2) is the volume of a d-dimensional sphere of unit diameter. In 
the second line of Eq. Q we have introduced the packing fraction = pvdcr'^ and the reduced virial coefficients 6„ = 
, Bn/{vd(j'^)"'~^ = 2('^"i)(""i)i3„/B2"^ The radius of convergence of the virial series (P), pconv = lim„^oo \BnlBn+i\, 
. is the modulus of the singularity of Z{p) closest to the origin in the complex p plane. If such a singularity were located 
' on the positive real axis, then all the virial coefficients i?„ would be positive for large n. 
5^ , The exact expression for the third virial coefficient is^'^ 



^3 ^ ^3/4(^/2 + 1,1/2) 
Bl S(d/2 + 1/2,1/2) 



where B{a,b) — r(a)r(6)/r(a + b) is the beta function and Bx{a,b) is the incomplete beta function. Both B2 
and i?3 are positive definite for arbitrary d. The analytic evaluation of the fourth virial coefficient is much more 
involved. Luban and Baram^ derived exact expressions for two of the three terms contributing to B4 and proposed 
a semi-empirical formula for the remaining contribution. More recently, Clisby and McCoy^Siffi showed that B4 can 
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be evaluated exactly for any even dimension d and gave the explicit results for d = 4, 6, 8, 10, and 12. They also 
computed numericalljJ^Siffi the fifth and sixth virial coefficients through d = 50. The results show that, while 
remains positive, B4 and become negative for d > 8 and d > 6, respectively. This suggests the possibility that, 
even in the three-dimensional case, there might be negative virial coefficients i3„ for sufficiently large nJ^- The fact that 
the virial coefficients are not positive definite and that they may have alternate signs is of importance in connection 
with the radius of convergence of the virial expansion as mentioned before. 



A. The scenario for high dimensionahty 

The high-dimensionality limit of hard hyperspheres has been the subject of several studies^iiSiSi. By means of 
asymptotic methods and heuristic arguments, Frisch and Percus^ were led to the following scenario in that limit: 

• The fourth virial coefficient is negative. Beyond that term, the virial expansion is an alternating series. 

• The virial series is convergent for p < 1, where p = 27/^/"^ is the scaled density per dimension. In terms of the 
packing fraction, the virial series converges for 77 < 77C01W = 2^'^. 

• In the density range p < I the second virial term dominates over the remaining ones, so that 



1 + ^2/9 = 1 + y. 



(3) 



• Even though the virial expansion does not converge for p > 1 (oscillatory divergence), the truncated series © 
remains valid within the interval 1 < p < (1 — e)/5o, where e — 0{d^^) and 
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• At the density p — pa an infinite compressibility spinodal appears, thus indicating a first-order transition to the 
high-dimensional solid. 

In an independent paper, Parisi and Slanina^ reached similar conclusions from a toy model based on simplified HNC 



equations. They obtained that, while in the limit d 
1 < p < Po = ^e/2, one has 



00 the EOS for p < 1 is given by Eq. 0, in the interval 
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where 



A(p) = 



2p2 
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k(p) being the solution t( 
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In^ =ln(l + Vl-K^ 
e 
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Note that, since Iuk < In (l + \/l — k^) — Vl — for < k < 1, one has limd^oo A(p) = for p < po. Although, 
strictly speaking, Eq. Q cannot be extended to p > po, Eq. @ suggests that lim^j^oo A(p) = cx) in that density 
domain, in agreement with the phase transition noted by Frisch and Percus^ 



B. Approximate equations of state 

As in two and three dimensions, one can make use of approximate schemes to represent the EOS of hard hyper- 
spheres. Several proposals have been made in the literature for the EOS based on the knowledge of the first few virial 
coefficients i2iiS4iiii24i^ For illustration, we review here a few of them making use of the first three virial coefficients. 
The extension to higher virial coefficients is straightforward. 
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1. Truncated virial series 

The first obvious choice is the truncated virial expansion 

^p.ojW = 1 + M + W- (8) 

As discussed above, this simple approximation becomes more and more accurate in the stable fluid domain as the 
dimensionality increases. In a way analogous to Eq. ^ it is possible to define a truncated expansion Z^^fl] from the 
knowledge of the first n + 1 virial coefficients. 

2. Fade approximants 

One can also construct Pade approximants of the form Z^n^m] from the first n + m + 1 virial coefficients. For 
instance, 

Z[i,i][r]) = r 7 : (9 

02 - 037] 

Z[o,2] iv) = [1 - b^v + [bl - b^W] . (10) 

3. Colot-Baus approximation 

Colot and Baus^iS proposed (truncated) rescaled virial expansions, where the series expansion of (1 — rj)'^Z{ri), 
rather than that of Z{ri), is truncated. Let us denote by Z^^j(?7) the truncated rescaled virial expansion that makes 
use of the first n + 1 virial coefficients. For example. 



,BC , ^ _ 1 + (&2 - d)r] +]63_M + d{d - l)/2]r]^ 

(1-77)" 



'^[2fi]\ni M — ~ • y^^i 



The pole of order d at the (unphysical) packing fraction 77 = 1 is suggested by the scaled particle theory. 

4. Maeso-Solana-Amoros-Villar approximation 

Maeso et al^ combined the advantages of Pade approximants and rescaled expansions by proposing Pade approx- 
imants for (1 — r])'^Z{rj). A rescaled Pade approximant constructed from the first n + m + 1 virial coefficients will be 
denoted here as Z^^^{ri). Thus, 



MSAV/ N _ 1 b2-d+[d{d+l)/2 + b2{b2-d)-b3]7] 



b2~ d--[b3~b2d + d{d^l)/2]rj 



(12) 



By construction, Z^^f [i]) = Zf^^^^{i]). 



5. Song-Mason-Stratt approximation 

Using simple arguments. Song et al.^"'^^ proposed the following generalization to d dimensions of the celebrated 
Carnahan-Starling (CS) EOS for three-dimensional hard spheres:^^ 



(13) 
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6. Luban-Michels approximation 

On a different vein, Luban and Michelsi^ wrote the compressibility factor as 

7 f\ 1 + [b3/b2 - CMh/bs] ry 

1 - C{v)[b4/b3)v + [Civ) - 1] {bi/b2W 

The knowledge of the function ^(77) is equivalent to that of Z{r]). However, C,{rf) focuses on the high density behavior 
of the EOS, since Eq. (|14|l is consistent with the exact first four virial coefficients, regardless of the choice of C,{r}). 
The approximation Civ) = 1 is equivalent to assuming a Pade approximant Z\2s]{'n)- Instead, Luban and Michels 
observed that the computer simulation data for d = 2-5 favor a linear approximation Qirf) = a + hr], with coefficients 
obtained by a least-square fit to the simulation results for each dimensionality. 



7. Percus-Yevick theory 

It is noteworthy that the Percus-Yevick (PY) integral equation can be solved analytically in odd dimensions, as 
first pointed out by Freasier and Isbister^ and, independently, Leutheusseriii The latter concluded that, in general, 
the problem reduces to an algebraic equation of degree d — 3. Following his procedure, however, we find that for 
d > 9 this is not so (see the Appendix) and our calculations suggest that such degree should rather be 2(''~'^)/^ for 
d > 3. In any case, in five dimensions one has to deal with a quadratic equationi^ and explicit expressions for the 
virial and compressibility routes to the EOS can be obtained.^" A simple analysis of the solution for c? = 5, that 
as far as we know has not been carried out before, shows that the virial route incorrectly gives a negative value 
for Be: B^^'^ / = ~ -0.00286. The compressibihty route yields B^^-^B^ = ~ 0.00146, while 

the correct value is iJg/Sf — 0.00094.^^'^° Both routes consistently predict that Bs is negative, with subsequent 
coefficients alternating in sign. On the other hand, the virial route gives values for the magnitude of i?„ {n > 8) 
increasingly larger than the compressibility route: Bf^^'^ / B^^~^ « 0.66 + 0.77n. The alternating character of the 
virial series predicted by the PY equation for d = 5 is due to a branch singularity located on the negative real axis 
at rebranch = — (9 — 5-\/3)/6 ~ —0.0566243. The radius of convergence r?^onv — 0.0566243 of the PY solution for d = 5 
is larger than the value 2~^ = 0.03125 extrapolated from the radius limd^oo Vcom = 2~'^, but is close to the estimate 
r/conv — 0.052 made by Clisby and McCoy on the basis of Monte Carlo evaluation of sets of Ree-Hoover diagrams 
All these estimates are sensibly smaller than the packing fraction rjf = 0.19 at which freezing occurs for d ~ 5.'^iii2& 



8. Generalized Camahan-Starling approximation 

As is well known, the CS EOS for three-dimensional hard spheres can be interpreted as a weighted average between 
the PY virial and compressibility routes: 

Zcsiv) = aZpY-dv) + (1 - a)ZpY-v{r]), (15) 

where a = |. Given that the PY equation can be solved for odd dimensions, it is then natural to speculate about 
whether or not the prescription H15|) . with an adequate choice of the mixing parameter a, keeps being reliable for d > 3, 
even though the internal inconsistency between both routes seems to increase dramatically with the dimensionalityi 
In the five-dimensional case, one of us^^ showed that the choice a — ^ leads to values of Zqs in excellent agreement 
with computer simulations.^ This suggested that the choice a = {d + l)/2d might provide a good description for 
d > 3. Note that, while Eqs. (|13|) and (|15|l coincide at d = 3, they differ for d > 3, so they generalize the original 
CS EOS along different directions. An alternative generalization of the CS EOS was made by Gonzalez et al>i& They 
proposed a simple ansatz for the direct correlation function c(r), which reduced exactly to the PY theory for d = 1 
and d = 3 and gave results very close to the PY theory for other dimensions. Their generalized CS EOS consisted of 
a weighted average between the virial and compressibility routes obtained from their theory with a mixing parameter 
a = 2(2d- l)/5d. 



C. Aim of the paper 



The aim of this paper is threefold. First, we present the explicit solution to the PY equation in the case of a 
seven-dimensional hard-sphere fluid following the procedure introduced by Leutheusseri4 This allows us to derive the 
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TABLE I: Values of Bn/ Bl^ ^ for n — 3-8, according to the virial route of the PY approximation, the compressibihty route 
of the PY approximation, the CS-like approximation HSU with a = |, and the known exact resultsi— 



3 0.2822265625 0.2822265625 0.2822265625 0.2822265625 

4 -7.499694824 x 10"^ 2.155081431 x 10"^ 1.670906279 x 10"^ 9.873(4) x 10"^ 

5 1.235022893 x 10"^ 5.116807918 x 10"^ 6.322378086 x 10"^ 7.071(7) x 10^^ 

6 -8.177005666 x 10"^ -1.865328120 x 10"^ -2.917274378 x 10"^ -3.52(2) x 10"^ 

7 6.553131160 x 10"^ 1.384246670 x lO'^^ 2.245727418 x 10"^ 

8 -5.762797816 x 10"^ -1.078783146 x 10"^ -1.859452258 x 10"=* 



EOS of the fluid both through the virial and the compressibility routes, as well as to analyze the behavior of the virial 
coefficients stemming out of them. As we will see, the singularity closest to the origin is again a branch point on 
the negative real axis, so the radius of convergence of the PY virial series is jyconv — 0.0100625. We conjecture that 
this value might be close to the (unknown) exact radius. Moreover, a Carnahan-Starling-like equation of state of the 
form 1)1511 with a = | is proposed. Secondly, we provide molecular dynamics results for the compressibility factor. 
To the best of our knowledge, this is the first time that simulation results are presented for hard hyperspheres in 
seven dimensions. The twenty densities considered range from the dilute regime (pu^ = 0.1 or r; = 0.0037) to our 
estimated freezing point {pa'' ~ 0.95 or 77 ~ 0.072). Finally, we perform a comparison between different proposals for 
the EOS of a seven-dimensional hard-sphere fluid with the simulation data. We observe that the proposals (|11() and 
(|13|l (which do not have any empirical parameter), 1)14)1 (which contains two fitting parameters), and 1)15)1 (with one 
fitting parameter) reproduce fairly well the simulation data. 

The paper is organized as follows. In Section^]we provide the solution of the PY equation for a seven-dimensional 
hard-sphere fiuid as well as the analysis of the virial coefficients arising from the derivation of the EOS using the virial 
and the compressibility routes. This is followed in Section Hill bv a description of the molecular dynamics simulation 
that was carried out to obtain the compressibility factor of the fluid. The results of the simulation are then used 
to assess the merits of various proposals that have been made in the literature for the EOS. The paper is closed in 
Section ll VI with further discussion of the results and some concluding remarks. 



II. SOLUTION OF THE PERCUS-YEVICK EQUATION FOR A SEVEN-DIMENSIONAL 

HARD-SPHERE FLUID 

As mentioned in SecU the solution to the PY equation for hard hyperspheres with d = odd reduces to an algebraic 
equation of degree 2^'^"^^/'^. The case d — b, which yields a quadratic equation, has been analyzed by several 
authorsi^ii^i^SiSi The highest dimensionality for which the algebraic problem certainly lends itself to an analytic 
solution is d = 7. A sketch of the general solution and some details of the particular cases d = 7 and d = 9 
are provided in the Appendix. It is shown there that the solution of the PY equation for seven-dimensional hard 
hyperspheres is given by the physical solution to the quartic equation l)A19p . In the Appendix it is also shown that 
for d = 9 the resulting algebraic equation is of eighth degree. 

A study of the solutions of Eq. IIA19)l shows that in the interval 0.446469 < ?/ < 1 the four roots are real. On the 
other hand, for < 77 < 0.446469 two of the roots become complex conjugates and only the other two roots remain real, 
the physical one being finite in the limit 77 ^ 0. The explicit solution to Eq. ()A19)) involves the term [P4(ri)PQ{r])]^^^ , 
where P^iTi) = 1 + 9^7^ + 202tj^ + ^r]^ + bOrj'^ and Peirj) ^ 1 + 99r/ - 2|Z^2 _ 339^3 _ m2^^4 _^ 695^5 ^ ^rj^. 
As a consequence, the solution possesses branch points at the zeroes of P4{ri) and Pe{ri). The zero of P4{ri) closest 
to the origin is il'hranch — —0.0108868, while that of P(){r]) is ?7branch — —0.0100625. Therefore, the radius of 
convergence of the virial series for a seven-dimensional hard-sphere fluid described by the PY approximation is 

Vconv = ^branch I ^ 0.0100625. 

Table HI gives the first few values of the PY virial coefficients obtained from the virial and the compressibility 
routes. As far as we know, the exact values B'^ of the virial coefficients of seven-dimensional hard spheres are 
known up to n = 6 onlyi ^^i^° They are listed in Table as well, which also gives the CS-like values B^^/B^-\ where 
B^s ^ aB^^-'' + {1- a)B^^-^ with the simple choice a = |. Note that the choice a ~ 0.6 would make ^ ~ B"^, 
whereas the choice a ~ 0.7 would make B'^^ ~ and Bq'^ ~ Bq^. However, comparison with molecular dynamics 
simulations (see Section 1111}) favors a ~ 0.8. 

From Tabled we observe that the virial route of the PY approximation incorrectly yields a negative value for the 
fourth virial coefficient (which actually becomes negative for d > g^SSiSSiSS) while the compressibility route predicts 
the correct sign.^^ We have computed B^^~^ and Bf^^~^ for values of n much larger than those displayed in Tabled 
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FIG. 1: Plot of the ratios |6„ /6„+i''| (diamonds), |b„ ""/fe^+l"! (triangles), and |fe„ /6„+i| (circles). The horizontal lines 
correspond to the apparent radius of convergence r^conv — 0.0108868 and to the true radius of convergence r;fonv — 0.0100625. 



The results indicate that sgn{B^^-^) = for 5 < n < 97 but sgn (bP^-v^ ^ f^j. ^ > gg. analogously, 

sgn(BPY-c) ^ for 5 < n < 80 but sgn = (-1)" for n > 81. Therefore, both routes synchronize 

their signs for 5 < rt < 80 and again for n > 98. This peculiar behavior of the alternating character of the 
virial series seems to be a consequence of the proximity between the two branch point singularities closest to the 
origin, Tybranch — —0.0100625 and ?7(,ranch — —0.0108868, both located on the negative real axis. To confirm this 
interpretation, we plot in Fig. [H the ratios \b^^~'' /b^J'^], Ibn^'"" /tH.'^], and Recall that the radius of 

convergence of the virial series is ?7conv = hm„^oo \bn/bn+i\- FigureQshows that for n < 50 the ratio |6„/5„+i| seems 
to converge from above to the apparent radius of convergence r?conv = "^branch — 0.0108868. However, the true radius 
^conv = — '7branch — 0.0100625 is reached from below for n > 100. 

As mentioned in Sec. the radius of convergence predicted by the PY approximation in the five-dimensional 
case is rj^^^ — 0.0566243. When going to the next odd dimensionality, the radius of convergence has shrunk to 
r]^^^ ~ 0.0100625. In terms of the scaled density per dimension introduced by Frisch and Percus,^^ the radius of 
convergence is p^o^v — 1-126 for d = 5 and Pconv — 1-037 for d = 7. Therefore, it nicely tends to converge to the 
expected value pconv = 1 as d oo. In addition, the PY value of r^conv = 2^~''i32 hm„^oo \Bn/ Bn+i\ ior d — 7 is 
consistent with estimates obtained from Table XVIII of Ref. By assuming that the Ree-Hoover ring diagrams 
dominate for high one has r/conv < S^^SalBg/Siol « 2-'^0.0132/0.0143 ~ 0.014 for d = 7, which agrees with 
the PY value ry^onv — 0.0100625. Chsby and McCoy's estimateSS ?7conv — 0.052 for d = 5 is also close to the PY 
value 77^onv — 0.0566243. All of this leads us to conjecture that the PY solution gives a fair estimate of the radius 
of convergence of the true virial series for high dimensionalities. Pushing this conjecture even further, we can expect 
the true radius of convergence to be due to a singularity (pole or branch point) located on the negative real axis, so 
that the virial coefhcients alternate in sign beyond a certain order. Figure 12 shows the virial coefhcients S^^'^/Bj 
B^^~'^ / B^ ^^ , and B^^/B2~^ in the seven-dimensional case. In the spirit of the above conjecture, one may speculate 
that the exact values of Bn/ Bl^^^ lie in between i?^^"^/i?2 ^^ and i?,^^"'^/i?2 perhaps not far from the interpolated 
values B'^^ / B2^^ . The reduced virial coefficients Bn/B2^^ start decreasing in magnitude, reach a minimum around 
n = 10, and then grow with n. The fact that the PY solution in the three-dimensional case does not possess a branch 
point singularity, so that all the virial coefficients remain positive, casts some doubts as to whether the true virial 
series fails to converge for densities close to the freezing density — 0.494. In any case, the true radius of convergence 
for d = 3 cannot be larger than the crystalline close-packing value r/cp — 7r\/2/6 ~ 0.7405, while the PY solution has 
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FIG. 2: Plot of the virial coefficients B'^''-yB^-' (diamonds), B^^'^/Bl^-' (triangles), and B^^/B^-' (circles). 



III. MOLECULAR DYNAMICS SIMULATIONS 
A. Method 



The numerical simulation was implemented by using the same algorithm as described in Ref. 25, which is also based 
on the work of Michels and TrappeniersSi and Luban and Michelsi^ for four- and five-dimensional hyperspheres. We 
are not aware of any previous computer simulation of hard hyperspheres of a dimension higher than d = 5. We have 
chosen the molecular dynamics method instead of the Monte Carlo method because that gives us the possibility of 
testing our code by applying it to d = 4 and d = 5 and comparing with the results of Refs. ISlllSl 

For our simulations, in order to keep the computing time within reasonable limits and at the same time being able to 
examine a wide density range, the initial configuration is chosen to be the one obtained by placing iV = 64 hyperspheres 
in a unitary cell of a d-type lattice. The simulation cell is a hypercube of side L and volume V = L''' — N / p and the 
minimum image convention and periodic boundary conditions in all directions have been applied, in the same way as 
in the 3D c&seM> 

During the simulations only binary collisions are taken into account, while collisions between three or four particles 
are ignored. The collision time for every pair of particles is calculated and the the smallest value is obtained. All 
the particles are moved during this time at constant velocity. The pair of particles that suffers a collision is treated 
according to impulsive dynamics and the velocities are changed; in this step the hard coUisional virial is calculated. 
This allows one to evaluate the excess compressibility factor as 



Z - 1 



1 



(16) 



where (w^) is the mean square velocity. At is the simulated time, Vij is the relative position vector between colliding 
particles i and j, and Av^ the change in velocity of the particle i on collision. 

The equation of state is achieved by changing the diameter a of the particles, in such a way that the reduced 
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density p* = pa changes, and letting the system to relax up to an equilibrium pressure. The errors associated with 
our calculation were computed following standard methods for errors in equilibrium averagesj21 

As mentioned above, before running the program for d = 7, it has been previously validated for d = 4 and d = 5, 
reproducing the excess compressibility factor obtained by Michels, Trappeniers, and Luban^*i^ 

B. Results 

We have computed the compressibility factor for densities 0.1 < p* < 1.90 with a step Ap* = 0.1, as well as for 
p* = 1.95. The simulation data obtained by our molecular dynamics simulations are listed in Table HTl At the largest 
density p* = 1.95 [r] = 0.0720) the compressibility factor presents a dramatic drop. We interpret this as an indication 
of the freezing transition. Consequently, the density at which the seven-dimensional fluid of hard spheres freezes 
can be estimated as p*j < 1.95 or, equivalently, 77/ < 0.072. From Fig. 5 of Ref . l2^ one can observe that ln?7/(d) is 
almost a linear function of the dimensionality d, with a slight negative curvature. According to this, knowing the 
freezing densities and r]f{d + 2), one can estimate the freezing density r]f{d + 4:) as ?7/(d-|-4) < rf^{d + 2)/r]f{d). 

Given that 77/(8) ~ 0.494 and ??/(5) ~ 0.19, one has 77/(7) < 0.19^/0.494 ~ 0.073, in close agreement with our 
estimate. An independent estimate based on a conjecture by Colot and Baus^ confirms again this value. These 
authors suggested that the ratio of length scales [?7/(d)/??cp('i)]^^'^ is practically independent of d, so that r]f{d+ 2) ~ 
Vcpid + 2)['rif{d)/r]cp{d)]'^'^'^^^^'^. The general expression for the close-packing fraction r]cp{d) is not known, but for 
d < 25 the values are not far from Blichfeldt's upper estimate^^ r]cp{d) < 2~''/^(d -|- 2)/2. Using ?7/(5) ~ 0.19, we get 
77/(7) ~ 0.076, which is again consistent with our estimate. 



TABLE IL Compressibility factor as a function of ry from the simulation data and for different approximations. The numbers 
in parentheses indicate the statistical error in the last significant digit. 
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3.55469 


3.51502 


3.55399 


3.55394 


0.0295 


4.0989(7) 


4.06234 


3.92342 


4.09012 


4.05480 


4.05948 


4.01769 


4.09372 


4.04037 


4.09264 


4.10109 


0.0332 


4.7013(5) 


4.64302 


4.44519 


4.68258 


4.64133 


4.65183 


4.57911 


4.68483 


4.61713 


4.68326 


4.70519 


0.0369 


5.389(1) 


5.27757 


5.00555 


5.33198 


5.28785 


5.30924 


5.18944 


5.33130 


5.24971 


5.32910 


5.36929 


0.0406 


6.051(1) 


5.96955 


5.60592 


6.04228 


6.00015 


6.04073 


5.85164 


6.03659 


5.94307 


6.03358 


6.09511 


0.0443 


6.8179(6) 


6.72276 


6.24777 


6.81775 


6.78456 


6.85731 


6.56896 


6.80433 


6.70273 


6.80033 


6.88228 


0.0480 


7.6325(1) 


7.54121 


6.93269 


7.66292 


7.64794 


7.77255 


7.34490 


7.63837 


7.53489 


7.63317 


7.72727 


0.0517 


8.5133(2) 


8.42922 


7.66233 


8.58259 


8.59769 


8.80333 


8.18328 


8.54278 


8.44645 


8.53613 


8.62220 


0.0554 


9.4294(3) 


9.39134 


8.43841 


9.58192 


9.64172 


9.97083 


9.08829 


9.52186 


9.44519 


9.51348 


9.55355 


0.0591 


10.492(1) 


10.4324 


9.26275 


10.6664 


10.7885 


11.3020 


10.0645 


10.5801 


10.5398 


10.5697 


10.5009 


0.0628 


11.570(3) 


11.5577 


10.1372 


11.8417 


12.0469 


12.8317 


11.1168 


11.7224 


11.7400 


11.7096 


11.4365 


0.0664 


12.694(1) 


12.7725 


11.0638 


13.1142 


13.4266 


14.6056 


12.2507 


12.9537 


13.0569 


12.9381 


12.3252 


0.0701 


13.907(3) 


14.0828 


12.0446 


14.4904 


14.9374 


16.6847 


13.4722 


14.2794 


14.5029 


14.2607 


13.1265 


0.0720 


9.03944(6) 


14.7756 


12.5559 


15.2196 


15.7454 


17.8637 


14.1178 


14.9794 


15.2786 


14.9589 


13.4810 



TableEJalso gives some theoretical values: the PY predictions Zpy-v and Zpy-c, the CS-like interpolation (|15|l with 
a = |, the truncated virial expansion .^[4,0] i the Fade approximants ^[2.2] and -^[3,2] [the three latter being obvious 
extensions of the approximations (|5|l- (|10|l ] , the rescaled virial expansion ^p^j defined by Eq. (|11|) , the rescaled Fade 
approximant Zj^^^^ defined by a natural extension of Eq. (|12() , the SMS approximation (|13|l , and the LM proposal 

Although the knowledge of the sixth virial coefficient Bq would allow one to consider the truncated series .^[5,0] j 
it is not included in Table ^ because it turns out to be clearly inferior to .^[4,0] ■ This is a consequence of the fact 
that Bq < 0, so that .^[5,0] < -^[4.0] j while for small and moderate densities .^[4,0] < -^simui- This is a strong indication 
that the unknown seventh virial coefficient B-^ must be positive. Among the different Fade approximants that can 
be constructed from the knowledge of the first six virial coefficients, the best agreement with the simulation data is 
presented by Z[2^2] for p* < 1.4 (77 < 0.0517) and by ^[3^2] for p* > 1.4 {rj > 0.0517) . It is interesting to note that 
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FIG. 3: Plot of the simulation values of the function (^{rj) defined by Eq. dJ. The dashed line is the linear fit 1)18^ . 



both Pade approximants have poles on the negative real axis (at r/ ~ —0.079 in the case of ^[2.2] and at 77 ~ —0.025 in 
the case of ^[3,2]), so that the extrapolated virial coefficients have alternating signs. Paradoxically, while the rescaled 
expansion .^p^j incorporates the first three virial coefhcients only, it exhibits a better agreement vifith simulation 
than those rescaled expansions that can be constructed with the first fom', five, or six virial coefficients, so the latter 
are not included in Table ITU Analogously, the best performance among the rescaled Pade approximants corresponds 
to ^plj^^. Interestingly, the SMS proposal [cf. Eq. and the approximation •^[^'^j yield practically equivalent 

results. The difference between both EOS is 



(1-77) 



y (35 - 35?/ + 21r;^ 



(17) 



This corresponds to a relative difference smaller than 0.14% for the density range considered in the simulations. 

Two of the theoretical EOS included in Table ITll namely Zcs a-nd Zlm, have an empirical character. The proposal 
(|15|1 is based on the observation that the two PY routes tend to bracket the simulation data, as happens in the three- 
dimensionalS and five-dimensional2i cases. We have found that the value a = | of the parameter is the simplest 
rational number that makes Zqs reproduce fairly well the simulation values. In the case of the Luban-Michels EOS 
(|14|1 one fits ({rj) to a linear function. Figure O shows the simulation values of C(^)- As in the five-dimensional 
case^ (^(77) is an increasing function of 77, while it is a decreasing function for d = 2-4. A linear fit in the interval 
0.5 < p* < 1.9 (0.0185 < ?7 < 0.0701) yields 



(^{j^) ^ -5.81 + 88.277. 



(18) 



The column labeled Zlm in Table UTI has been evaluated using the fit (|18|l . On the other hand, our simulation data 
in Fig. 13 seem to indicate a negative curvature of C('7)- 

Table shows that up to p* — 0.8 (77 = 0.0295) all the different theoretical results tabulated, including the 
simple truncated virial expansion .^[4.0], behave relatively well. For larger densities, .^[4.0] tends to overestimate the 
simulation data, while the Pade approximants .^^[2,2] and .^[3,2] tend to underestimate them. The best global agreement 
is presented by Zcs, -Z^lm, Z\2%]^ ^^'^ ^sms- This is especially noteworthy in the case of the two latter approximations, 
since they do not contain fitting parameters and, moreover, only the knowledge of the first three virial coefficients is 
exploited. This contrasts with Zlm, which includes the fourth virial coefficient and contains two fitting parameters. 
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FIG. 4: Compressibility factor as a function of the packing fraction. The circles represent simulation data, the dashed line 
represents Zpy-c, the dash-dotted line represents Zpy-v, and the solid line represents the CS-like interpolation H15|l with a — ^. 

On the other hand, Zqs belongs in a different class of approximations. Given the involved algebraic structure of 
the PY solution, Zcs does not intend to represent a practical recipe to the EOS of a seven-dimensional hard-sphere 
fluid. Instead, its role is to highlight the fact that the two PY routes keep bracketing the simulation data, so that an 
interpolation between them with a density-independent parameter a is rather accurate, as graphically illustrated in 
Fig. ^ This gives some confldence on the expectation that some of the analytical properties of the PY solution (e.g., 
alternating character of the virial series, branch points located on the negative real axis, . . . ) may shed light on the 
true behavior of the exact series. 



The results of the previous sections deserve further discussion. To begin with, to our knowledge this is the first time 
that a molecular dynamics simulation has been carried out on a seven-dimensional hard-sphere fluid. The simulation 
strategy that we adopted implied a compromise between computer process time and density range to be explored 
and the outcome is rather encouraging. The availability of simulation data for the EOS of the fluid allowed us to 
locate the freezing transition and also to assess the merits and limitations of various proposals that have been made 
in the literature for the compressibility factor of hard hyperspheres. From this analysis it is clear that even simple 
approximations such as Zj^^^^ and .^sms do a reasonably good job and that, as it occurs in other dimensionalities, the 
virial and compressibility routes to the EOS in the PY approximation keep bracketing the simulation data, so that 
a Carnahan-Starling-like recipe of the form of Eq. H15() turns out to be rather accurate. However the parameter a 
seems not to follow a simple relation as the ones suggested by Gonzalez et aliiSiii or SantosiSi 

We also presented the explicit solution of the PY equation for a hard-sphere fluid in 7D. Such a solution allowed 
us to carry out an analysis of the virial coefficients arising both in the virial and in the compressibility routes and 
to determine the radius of convergence of both virial series. The results indicate some peculiar behavior of the virial 
coefficients with the virial route incorrectly predicting a negative fourth virial coefficient. The radius of convergence 
of the virial series is due to a singularity (branch point) located on the negative real axis and therefore what one has 
is an alternating series. Because of the good agreement between our value of the radius of convergence of the virial 
series and other independent estimates and the similar results obtained for d = 5, it is tempting to conjecture that the 
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PY solution for even higher dimensionahties should provide a rather accurate estimate of the radius of convergence 
of the true virial series and that it is the existence of singularities on the negative real axis (either poles or branch 
points) which determines such radius. 

As a final point it is worth commenting that in this case our analysis was facilitated by the fact that we could 
combine both the analytical and the simulation results. And due to the common features such as the freezing transition 
that hard-core systems in different dimensionalities share, the expectation and the hope is that the present results 
shed some more light on the thermodynamic properties of such systems. As far as the high dimensionality limit is 
concerned, our results provide some support to the scenario of Frisch and Percus mentioned in the Introduction in 
the following sense. The solution to the PY equation predicts an alternating virial series. Further, the values of the 
scaled density p that one obtains for the radius of convergence (p = 1.13 for d = 5, p = 1.04 for d — 7, and the 
number p ~ 1.02 coming out of our preliminary calculations for d ~ 9) are consistent with a limiting value of p = 1 
for d ^ oo. Also, the fluid range in c? = 7 is reasonably well accounted for by the first three or four virial coefficients 
so that it is conceivable that for infinite dimensionality only the second virial coefficient will be the dominant term. 



Acknowledgments 

M.R. acknowledges the financial support of CONACyT through project No. I38644-E. M.L.H. acknowledges the 
hospitality of Universidad Complutense de Madrid (Spain), where the final version of the paper was prepared, as 
well as the financial support of DGAPA-UNAM during his sabbatical stay in Madrid. The research of A.S. has been 
partially supported by the Ministerio de Ciencia y Tecnologia (Spain) through grant No. BFM2001-0718 and by the 
European Community's Human Potential Programme under contract HPRN-CT-2002-00307, DYGLAGEMEM. 



APPENDIX A: SOLUTION OF THE PERCUS-YEVICK EQUATION FOR HARD HYPERSPHERES 

For simplicity, in the remainder of this Appendix we set a = 1. 

In the Percus- Yevick approximation, the structure factor S{q) of a hard-sphere fluid in d = 2fc + 1 dimensions is 

S{q) = ^ 1 , (Al) 

Q{q)Q{-q) 

where 

Q(g) = l-A/ dre*«''Q(r), X = {2tt)'' p = 2^''{2k + ly.lrj, (A2) 

Q{r) having the form 

Q(,)./ELoQ™('^-l^+^ o<r<i, (^3) 

[^0, r ^ 1. 

The k -\- 1 coefficients {Qm} a^i^c functions of the density determined by the two hnear equations 

fc 

(-1)^ -m'^Qk + A ^ (-1)- ^"^ , A;>0, (A4) 

m— 



(-1)^ = -(fc - m'-'Qk-i + A ^ (-1)- ^, ^ , fc>i, (A5) 

m— 

plus the k — 1 nonUnear equations 

m — 1 

-A^(-1)''Q('^Ho)Q(2™-'')(o), 0<m<fc-2. (A6) 

Here Q^'^^{r) represents the v-th derivative of the function Q{r). For fc = (d = 1), Eq. (jA4|) gives the exact solution 
for hard rods. For fc = 1 (d = 3), Eqs. (|A4|) and ljA5|l are sufficient to find the solution of the PY equation. However, 



g(2™+i)(0) = -A(-l)"+i Q(")(0) 
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for fc > 2 (d > 5) one needs in addition Eq. ljA6|l . so that the problem reduces to solving an algebraic equation which, 
as we will argue below, is likely to be of degree 2^~^ — 2^'^"^^/^. 
In the limit 7y ^ 0, it is easy to verify that 

\imQ„, = {-lf+'—( , limg(r) = (-l)'=+i — (r^-lf, (A7) 
»?^o fc! \mj ri^o k\ 

lim g(2™)(0) = („l)r»+l2-fc (^H! g(2m+l)(o) ^ q (^g) 

m!(fc-m)! j;^o 
In general, one can expand the coefficients Qm in powers of tj: 



Qm(^) (A9) 

where Qm.o is given by the first equation of l|A7p . Of course, the full nonlinear dependence of the coefficients Qmiv) 
can be obtained from the solution to the set of equations l)A4p - lA6p . either analytically (fc < 3) or numerically (/c > 4). 

Once one has determined the fmiction Q{r), the structural properties of the fluid are given by Eqs. IjAip and IIA2|l . 
In particular, the long wavelength limit of the structure factor and the contact value of the radial distribution function 
are, respectively, 

S{q = {)) = i 2, (AlO) 

5(1+) = (-l)'=+ifc!Qo. (All) 
^ = l + 2''-i77.g(l+), (A12) 



The virial route to the EOS is given by 



while the compressibility route is 



X = kBT[^\ ^S{q^O). (A13) 



dp 



T 



Inserting the expansion l|A9|l into Eqs. (|A12p and l|A13|l we get the virial coefficients along both routes: 

n 

€17 = 22'=(-l)'=+lfc!Qo,„, hll\^ - 22'=(fc!)2-_- ^ Qk,„^Qk.^-m. (A14) 



m— 



1. The case d — 7 



Now we particularize to the seven-dimensional case {k = 3), the unknowns being Qm, m = 0, 1,2,3. Since the 
two nonlinear equations (|A6|) involve the derivatives Q^™"* = Q^"^H^), it is more advantageous to work with the set 
rather than with the set {Qm}- The latter can be expressed in terms of the former as 



/Qo\ 


( 20 


10 


2 


'A 


( \ 


Qi 


45 


25 


11 

2 


! 




Q2 ~ 


36 


21 


5 






\qJ 




6 


3 
2 


11 





(A15) 



Equations HA4|I and HA5|I . plus Eq. HA6|1 with to = yield 



Qfi' = -33607yQ(°)', 



(A16) 
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TABLE III: 


Values of the coefficients Q 


,„^„, defined by Eq. IIA9II. for n = 0-6. 




n 


Qo.n 


Ql.n 


Q2,n 







0.166666667 


0.250000000 


0.125000000 


0.020833333 


1 


3.010416667 


7.888020833 


5.812500000 


1.333333333 


2 


-5.119791667 


-59.313802083 


-41.072916667 


-6.541666667 


3 


5.395897352 x 10^ 


4.247567790 x 10^ 


3.201932292 x 10^ 


6.540590278 x 10^ 


4 


-2.286456505 x 10* 


-2.488562475 x 10^ 


-1.884527380 x 10^ 


-3.841568446 x lO" 


5 


1.172728501 X lO'* 


1.635350292 x lO'^ 


1.241794603 x lO'^ 


2.538798289 x 10^ 


6 


-6.600274174 x 10^ 


-1.143524052 x lO'' 


-8.688923174 x 10* 


-1.778764760 x 10* 



Q 



(2) 



1 + 96Q(°) {1 - 5?7 [3 + 112Q(") (3 - lOry)] } 



1577 + 192Q("' |2 



8(1-77) 

53 + 280Q(°)(3 - 1077)^ - IOO77 I 



(A17) 



(A18) 



Thus, the parameters Q'^^\ Q'"^\ and Q*^'^^ are given as exphcit quadratic functions of Q'-*'-' — Q(0). Finally, insertion 
of Eqs. (|A16p - (|A18|l into Eq. ljA6|l with m = 1 leads to the quartic equation 



8 - 15r; + 192Q(°) |2 - 77 188 - 135r/ + 1960g(°) 
X (3 - 1077 (i-84Q(")(1-77)))]]}} =0. 



3 - 4?7 9 - 10?7 + 240Q(") (1 - r/) 



(A19) 



Although an explicit expression exists for the physical root of Eq. (|A19p , it is of course too cumbersome and will be 
omitted here.^^ 

Table Hill shows the first few coefficients Qm,n- The exact values are rational numbers, but they become more and 
more involved as the order n increases and so they are expressed in real form in Table ITTll From Eq. (|A14|) we can 
obtain the virial coefficients corresponding to the virial and the compressibility routes. The first few values are listed 
in Table n 



2. The case d = 9 



We will now sketch the result for the case d = 9 following the same procedure. For fc = 4, the set {Qm} can be 
expressed in terms of the set {Q*-™'} as 



/Qo\ 

Qi ^ 

Q2 

Q3 



( 70 


35 
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119 


27 
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2 
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90 


22 


35 


20 


5 



24 



e ! 



(A20) 



In addition, the fifth derivative Q^^^ is a linear combination of {Qm} and hence of the first four derivatives: 

g^'^' = -20 (336g(°^ + 2iog(^^ + eog^^^ + log^^^ + g^''^) . (A2i) 

The nonlinear equations l|A6|l with m = 0, 1, 2 allow one to express the odd derivatives in terms of the even ones as 

Q(i) ^ _^Q(o)2, (A22) 



g 



(3) 



A3 



1(0)^ 



(o)n(2) 



Ag^'"g 



(A23) 
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Q(^) - -T?Q*°^' + ^Q^"^'Q^'^ - ^Q^'^' - AQ(")g(^), (A24) 
Id 2 2 

where A = 241 920?7. Next, insertion of Eqs. HA22p and (|A23p into the hnear equations (|A4|) and (|A5|) yields Q*^^ and 
Q'-^^ as nonlinear functions of Q'"'. FinaUy, by equating the right-hand sides of Eqs. (jA21|) and (|A24|) one gets a 
closed algebraic equation of eighth degree for Q'^^K A preliminary analysis of this equation indicates that its physical 
solution possesses a branch point at ?7branch — —0.0023945, so that the radius of convergence of the PY virial series 

would be 77conv = |?7branch| — 0.0023945. 

We have checked that for d = 11 the resulting equation is of degree 16. Therefore, it seems plausible that in the 
general case d = 2fc + 1 the degree of the equation for Q^°^ is 2*^"^ = 2^'^~^^^^. 
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